function pde = fourCurl3data
%% FOURCURL3DATA data for fourth order equation involving curls
%
%   u = [(x^2-1)^2*(y^2-1)^2*(z^2-1)^2 
%              0
%              0];
%
%   curlu = [0
%               (4z^3-4z)*(x^2-1)^2*(y^2-1)^2
%               -(4y^3-4y)*(x^2-1)^2*(z^2-1)^2 ];
%
%   curlcurl u = [-(12y^2-4)*(x^2-1)^2*(z^2-1)^2-(12z^2-4)*(x^2-1)^2*(y^2-1)^2
%                      (4y^3-4y)*(4x^3-4x)*(z^2-1)^2
%                      (4z^3-4z)*(4x^3-4x)*(y^2-1)^2];
%
%   curlcurlcurlu = [0
%                          -(4z^3-4z)*(12x^2-4)(y^2-1)^2-(4z^3-4z)*(12y^2-4)*(x^2-1)^2-24z*(x^2-1)*(y^2-1)^2,
%                          +(4y^3-4z)*(12x^2-4)(z^2-1)^2+(4y^3-4z)*(12z^2-4)*(x^2-1)^2+24y*(x^2-1)*(z^2-1)^2]
%
%   curlcurlcurlcurl u=[(12y^2-4)*(12x^2-4)(z^2-1)^2+...
%                               24*(x^2-1)^2*(z^2-1)+...
%                              (12y^2-4)*(12z^2-4)*(x^2-1)^2+...
%                              (12z^2-4)*(12x^2-4)(y^2-1)^2+...
%                               24*(x^2-1)^2*(z^2-1)+...
%                              (12z^2-4)*(12y^2-4)*(x^2-1)^2,
%                               
%                             ?24x*(4y^3-4y)*(z^2-1)^2-...
%                               24y*(4x^3-4x)*(z^2-1)^2-...
%                              (4y^3-4y)*(12z^2-4)*(4x^3-4x),
%
%                              ?24x*(4z^3-4z)*(y^2-1)^2-...
%                               24z*(4x^3-4x)*(y^2-1)^2-...
%                              (4z^3-4z)*(12y^2-4)*(4x^3-4x)]
% 
%   f = curlcurlcurlcul u + u;
%
%  - u satisfies zero Dirichlet boundary condition at cube [-1,1]^3
%
% Added by Lin Zhong, June 2013.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

pde.f = @f;
pde.exactu = @exactu;
pde.curlu = @curlu;
pde.g_D = @exactu;
pde.curlcurlu = @curlcurlu;
pde.curlcurlcurlu = @curlcurlcurlu;
pde.g_N = [];

    function s = exactu(p)
        x = p(:,1); y = p(:,2); z = p(:,3);
        s = [(x.^2-1).^2.*(y.^2-1).^2.*(z.^2-1).^2 ,0*y , 0*z];
    end

    function s = curlu(p)
    x = p(:,1); y = p(:,2); z = p(:,3);
    s= [0*x ,     (4*z.^3-4*z).*(x.^2-1).^2.*(y.^2-1).^2 , -(4*y.^3-4*y).*(x.^2-1).^2.*(z.^2-1).^2 ];
    end

    function s = curlcurlcurlu(p)
    x = p(:,1); y = p(:,2); z = p(:,3);
    s = [0*x,  -(4*z.^3-4*z).*(12*x.^2-4).*(y.^2-1).^2-...
        (4*z.^3-4*z).*(12*y.^2-4).*(x.^2-1).^2-24*z.*(x.^2-1).^2.*(y.^2-1).^2,+(4*y.^3-4*y).*(12*x.^2-4).*(z.^2-1).^2+...
        (4*y.^3-4*y).*(12*z.^2-4).*(x.^2-1).^2+24*y.*(x.^2-1).^2.*(z.^2-1).^2];
    end

    function s = f(p)
        x = p(:,1); y = p(:,2); z = p(:,3);
        s=[(12*y.^2-4).*(12*x.^2-4).*(z.^2-1).^2+...
                              24*(x.^2-1).^2.*(z.^2-1).^2+...
                             (12*y.^2-4).*(12*z.^2-4).*(x.^2-1).^2+...
                             (12*z.^2-4).*(12*x.^2-4).*(y.^2-1).^2+...
                              24*(x.^2-1).^2.*(y.^2-1).^2+...
                             (12*z.^2-4).*(12*y.^2-4).*(x.^2-1).^2,     -24*x.*(4*y.^3-4*y).*(z.^2-1).^2-...
                              24*y.*(4*x.^3-4*x).*(z.^2-1).^2-...
                             (4*y.^3-4*y).*(12*z.^2-4).*(4*x.^3-4*x),  -24*x.*(4*z.^3-4*z).*(y.^2-1).^2-...
                              24*z.*(4*x.^3-4*x).*(y.^2-1).^2-...
                             (4*z.^3-4*z).*(12*y.^2-4).*(4*x.^3-4*x)];                    
        s = s + exactu(p);
    end

    function s = curlcurlu(p)
        x = p(:,1); y = p(:,2); z = p(:,3);
        s= [-(12*y.^2-4).*(x.^2-1).^2.*(z.^2-1).^2-...
            (12*z.^2-4).*(x.^2-1).^2.*(y.^2-1).^2, (4*y.^3-4*y).*(4*x.^3-4*x).*(z.^2-1).^2, (4*z.^3-4*z).*(4*x.^3-4*x).*(y.^2-1).^2];
    end
end